clear
format long



%--------------------------------------------------------------------------
%----------------------NOTES-----------------------------------------------
%--------------------------------------------------------------------------

%NAME: CRTBPpoincareGrid.m

%Dependencies:
%This program depends on 'CRTBP.m', 'BackwardCRTBP.m', 'magnitudeVelocity.m',
%'librationPoints.m', 'jacobiConst.m', 'interpCrossingData.m', and 
%'CRTBPpoincareNewton.m'.  

%All of these are of these are called by the present program except
%'BackwardCRTBP.m' which is called by 'CRTBPpoincareNewton.m'.

%USE: This program builds a Poincare map, where the grid of initial 
%conditions is specified by the user.

%Below are several USER DEFINED CONTROL PARAMETERS.  The user should
%set these (or leave them set to their default values) and save the
%program.  Then run in UNIX background mode.  If you run it from the MatLab
%command line as it is now, then when it finishes (a day or two later) it
%will not produce any output.  

%However you can type 'who' at the command prompt and see that all the
%variables are there.  To produce a Poincare section just plot the first
%two columns of 'moser' against oneanother.  Plot as dots instead of
%connecting with points.

%OUTPUT:
%The output a file containing the variable 'PoincareMap'.  This is the date
%that is usually plotted.  The proformancs variables are writen as well.

%PURPOSE: Similar to 'CRTBPpoincare.m' and it's decendants, except here
%instead of taking initial conditions all on the x-axis with velocity
%perpendicular to the axis, the program builds a grid, specified by the
%user.



%--------------------------------------------------------------------------
%--------------------USER DEFINED CONTROL PARAMETERS-----------------------
%--------------------------------------------------------------------------

%The parameter 'G' is available in case non-normalized units are
%desired.  When this is not the case just leave 'G' set to one.
        G=1;

%The value of 'mu' depends on the mass ratio fo the primaries.
%If mu is zero the CRTBP reduces to the Kepler problem in non-inertial
%(rotating) coordinates.  The coordinates are normailized so that mu is
%always less than or equal to 1/2, with 1/2 the case of equal masses.
%Other important values are
%
%mu=3.054248395e-06      (for the sun/earth system)
%mu=0.012277471          (for the earth/moon system)
%
%The value of mu must be set by the user (default is earth/moon).
    mu=0.012277471

%Once 'G' and 'mu' are defined we can compute the location of the libration
%points and their Jacobi Integrals.  These will be handy when we set the
%value of the Jacobi Integral for the simulation below.

    %AXULARY CONSTANTS:
    %Compute the location of the five libration points as points in three
    %space.  These can be used to set the value of the Jacobi Constant.
    [L1, L2, L3, L4, L5]=librationPoints(mu);

    %Output location to the screen;
    outL1=L1;
    outL2=L2;
    outL3=L3;
    outL4=L4;
    outL5=L5;
    %Now compute the Jacobi Energy at each collinear libration point

        %The velocity at any equilibrium is zero.
        v_equil=[0; 0; 0];

    %Compute the energy at each libration point
    CL3=jacobiConst(L3, v_equil, mu);
    CL1=jacobiConst(L1, v_equil, mu);
    CL2=jacobiConst(L2, v_equil, mu);
    CL4=jacobiConst(L4, v_equil, mu);
    %CL5=CL4 so no extra computation is needed

%The value of 'C' is the value of the Jacobi Integral that you want to
%compute a Poincare section for.  For values less than CL1 but greater than
%CL2 the 'neck' at L1 is open but the 'neck' at L2 is closed.  Then we can
%have transport between the earth and the moon, but all orbits (in the
%Hill's Region of the Earth and Moon) are bounded.  Then, while the user
%can set 'C' to anything they like, our default value of 'C' is:

    C=CL1-0.5*(CL1-CL2)    %0.5 here is arbitrary.  Any coefficent between
                           %zero and one (inclusive) will have the property
                           %described above.

    %A secondary default value is provided below.  At this energy the neck at
    %L1 is not open.  Hence for initial conditions begining inside the Hills
    %region of the primary, the time scale problems are avoided

    %C=1.001*CL1;   %(Comment out this line if using primary default)
                           

%The variable 'k' determins how many initial conditions will be used to
%generate the x axis of the Poincare map.
    k=300

%The variable 'l' determines how many initial conditions will be used to 
%generate the y axis of the Poincare map.
   l=1
    
%The variable 'itterates' decides how many crossings of the xz-plane will
%be computed for every initial condition.  This is equivalent to specifying
%how many itterates of the poincare map to compute for every point. Values
%near 1000 are best, but involve long runtimes.
    iterates=1000


%The variable 'x_begin' sets the left endpoint for the interval of points
%that will be itterated.  Want to stay away from -mu and 1-mu as these are
%singularities.
    x_begin=0.085

%The variable 'x_end' is the right end-point.  The determines how far to
%go to the right of 'x_begin'.
    x_end=0.95

%Determined the lower limit of the 'y grid' (which is really an xdot axis).
   xdot_lower=0.0
   
%Determin the upper limit of the 'y grid' 
   xdot_upper=0.5;
    
    

    
    
%Secondary Controls:
%NOTE: The following parameters do not need to be changed (they were fine
%tuned in the developement of the program.  Changing them should not hurt
%and could even help, but is not necessary).

    %sets the tolerence for the map.  The poincare map looks for crossings
    %of the xy plane.  epsilon determines how close a trajectory point must
    %be to the plane to be considered "on" the plan.  Standard value is
    %ten to the minus six.
        epsilon=1e-6

%The program integrates for a time, then checks for crossings.  It repeats 
%process untill it has found the desired number of crossings.  The variable
%below defines how long to integrate before looking for crossings.  Thsi is
%somewhat arbitrary.

   timeStep=10;
        
        
        
%--------------------------------------------------------------------------



%NOTE: The following parameters are necessary but depend on the above

    %if k=1 no steps are necessare and the following variable makes
    %no sense.
    if k>1
        %determine the step size between successive initial conditons
        xStep=(x_end-x_begin)/(k-1);
    else
        %if k=1 just set step to zero so it is at least defined later
        xStep=0;
    end %end seeing if steps are necessary

   %if l=1 no vertical steps are necessary
   if l>1
       xdotStep=(xdot_upper-xdot_lower)/(l-1);
   else
       xdotStep=0;
   end %end of xdotStep    
    
    
    
%INITIALIZE:
    %This is the number of iterates of the current initial conditions
    %found this far through the while loop
    x_nIterates=0;
    %This is the total number of points in the map so far
    totalIterates=0;
    %Local name of initial conditions:
    initial_n=0;
    %number of times through the while loop
    checkWhile=0;
    %number o times the intorpolator is called
    interpTimes=0;
    ReIntegrate=0;
    checkTotal=0;
    calledNewton=0;

%--------------------------------------------------------------------------
%----------------------COMPUTATION OF POINCARE MAP-------------------------
%--------------------------------------------------------------------------



for m=1:l
%this loop sets the vertical grid

for n=1:k
    %I like to see what loop I'm in at command line, just to see how fast
    %the program is running and where it is.  Comment this out or delete it
    %if you prefer.  Same goes for any 'check' variable
    check_n=n
    checkTotal=totalIterates
    %First the position
            x0=x_begin+(n-1)*xStep;  %Step to the next initial x0
            y0=1e-12;  %Make close to zero but with positive sign
            z0=0;
        %Then the velocity
            %Direction of Velociety
            v_x0=xdot_lower+(m-1)*xdotStep;
            v_y0=1.0;   %Take ezch initial velocity  normal to x axis
            v_z0=0.0;
            v0=[v_x0; v_y0; v_z0];
            %make it a unit vector
            u_v=v0/norm(v0);
            %the magnitude of v is given by the Jacobi integral
            V=magnitudeVelocity([x0; y0; z0],C, mu);
            %then the initial velociety is the unit vector times the magnitude
            xdot0=V*u_v(1,1);
            ydot0=V*u_v(2,1);
            zdot0=V*u_v(3,1);
    %Set the initial condition at the begining of the nth time
    %through the for loop (nth initial condition) using the data above:
            initial_n=[x0;
                       y0;
                       z0;
                       xdot0;
                       ydot0;
                       zdot0]';

    %THE MAIN WHILE LOOP:  This computes the number of times the
    %trajectory with initial condition 'initial_n' crosses the poincare
    %section with positive y velocity. The desired number of such crossings
    %is given by the integer 'iterates'.  The count of how many have been
    %found so far is in the raviable:
              x_nIterates=0;   %This should be set back to zero before
                               %every run through the while loop.

    %Begining of while loop
    while x_nIterates <=iterates
        checkWhile=checkWhile+1;
               
        %Intigrate this for time 'timeLength'.
           tspan=[0 timeStep];
           %numerical integration
           options=odeset('RelTol',1e-13,'AbsTol',1e-22);     %set tolerences
           [t, trajectory_n]=ode113('CRTBP',tspan,initial_n,options,[],G,mu);

       temp=size(trajectory_n);
       numSteps=temp(1,1);
           
           
       %Now check for pairs of points along the trajectory between which
       %the sign of y changes and the sign of ydot is positive
       for i=2:numSteps
           %Check for sign change of y and correct sign of ydot.
           if (sign(trajectory_n(i,2)) ~= sign(trajectory_n(i-1,2))) & (trajectory_n(i,5) > 0) 
                %It could happen that one of these two points is already a
                %crossing to the desired tolerence.  We check for this:
                if (abs(trajectory_n(i-1,2))<epsilon)
                    intersection=trajectory_n(i-1,:);
                elseif (abs(trajectory_n(i,2))<epsilon)
                    intersection=trajectory_n(i,:);
                else
                %If (as is most likely the case) neither of the points 
                %for which the sign change occurs are in tolerence
                %to be considered as crossings and we are farther than 'delta' 
                %from a singularitys, then a Newton method determines the 
                %crossing to the desired accuracy.
                    intersection=CRTBPpoincareNewton(t(i-1),trajectory_n(i-1,:),t(i),trajectory_n(i,:),epsilon,G,mu);
                    calledNewton=calledNewton+1;  %(Print to scren)
                end %end of the intersection finding conditional.

                %Keep track of how many iterates of the current initial
                %condition have been found
                    x_nIterates=x_nIterates+1;  %(print to screne)
                %Keep track of how many total intersections have been found
                    totalIterates=totalIterates+1; %(print to screne)
                %Stor the intersection information for plotting later
                    PoincareMap(totalIterates, 1)=intersection(1);
                    PoincareMap(totalIterates, 2)=intersection(4);
           end %end of pair checking if statement
       end %end of pair checking loop (i=2:numSteps index)

        %make sure variables are ready for next while loop: unless all the
        %iterates have been found we will go fhrought the while loop again.
        %In that case the new initial condition is the final condition of
        %trajectory_n

        initial_n=trajectory_n(numSteps,:);

    end %end of main while loop (x_nIterates <= iterates)

    %The set up for the next time through the 'n' loop is done at the
    %begining of the loop instead of here at the end.
    %So just end and go back.


end %end of n=1:k loop

end  %end of outmost for loop m=1:l 

PoincareGrid=PoincareMap

save closeData2 PoincareGrid



